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(57) ABSTRACT 

This invention presents a method for computing Instanta- 
neous Frequency by applying Empirical Mode Decomposi- 
tion to a signal and using Generalized Zero-Crossing (GZC) 
and Extrema Sifting. The GZC approach is the most direct, 
local, and also the most accurate in the mean. Furthermore, 
this approach will also give a statistical measure of the 
scattering of the frequency value. For most practical appli- 
cations, this mean frequency localized down to quarter of a 
wave period is already a well-accepted result. As this 
method physically measures the period, or part of it, the 
values obtained can serve as the best local mean over the 
period to which it applies. Through Extrema Sifting, instead 
of the cubic spline fitting, this invention constructs the upper 
envelope and the lower envelope by connecting local 
maxima points and local minima points of the signal with 
straight lines, respectively, when extracting a collection of 
Intrinsic Mode Functions (IMFs) from a signal under con- 
sideration. 

14 Claims, 7 Drawing Sheets 
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Figure 1(b) 
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COMPUTING FREQUENCY BY USING 
GENERALIZED ZERO-CROSSING APPLIED 
TO INTRINSIC MODE FUNCTIONS 

ORIGIN OF INVENTION 5 

The inventor of the invention described herein is an 
employee of the United States Government. Therefore, the 
invention may be manufactured and used by or for the 
Government for governmental purposes without the pay- 10 
ment of any royalties thereon or therefor. 

BACKGROUND OF THE INVENTION 
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vice-versa) as it crosses the zero axis. In the field of signal 
processing, for a subject signal, the zero crossing points are 
detected and counted for analyzing the signal or making a 
decision for further actions. For example, zero crossing 
points are used to define bandwidth of a signal, which is 
assumed to be stationary and Gaussian. The bandwidth can 
be defined in terms of spectral moments as follows. The 
expected number of zero crossings per unit time is given by 



The need for instantaneous frequency (IF) is a real one for 15 
data from non-stationary and nonlinear processes. If the 
process is non-stationary, the frequency should be ever 
changing, albeit at a slow rate. Then, there is a need for 
frequency value as a function of time, for the value will not 
be constant throughout. For the nonlinear cases, the fre- 20 
quency is definitely modulating not only among different 
oscillations, but also within one period. A detailed explana- 
tion of IF is disclosed in related patent application Ser. No. 
10/615,365, filed on Jul. 8, 2003, entitled “Computing 
Instantaneous Frequency by Normalizing Hilbert Trans- 25 
form”, inventor Norden E. Huang, which is incorporated by 
reference and assigned to the same assignee as this appli- 
cation. A computer implemented method of computing IF is 
disclosed in U.S. Pat. Nos. 5,983,162, 6,381,559, 6,311,130, 
all of which are also incorporated by reference. 30 

The above -disclosed method of computing IF includes 
two essential steps and the associated presentation tech- 
niques of the results. The first step is a computer imple- 
mented Empirical Mode Decomposition to extract a collec- 
tion of Intrinsic Mode Functions (IMF) from nonlinear, 35 
nonstationary signals. The decomposition is based on the 
direct extraction of the energy associated with various 
intrinsic time scales in the signal. Expressed in the IMF’s, 
they have well-behaved Hilbert Transforms from which 
instantaneous frequencies can be calculated. The second 40 
step is the Hilbert Transform of the IMF. The final result is 
the Hilbert Spectrum. Thus, the method can localize any 
event on the time as well as the frequency axis. The 
decomposition can also be viewed as an expansion of the 
data in terms of the IMF’s. Then, these IMF’s, based on and 45 
derived from the data, can serve as the basis of that expan- 
sion. The local energy and the instantaneous frequency 
derived from the IMF’s through the Hilbert transform give 
a full energy-frequency-time distribution of the data, which 
is designated as the Hilbert Spectrum. 50 

However, there is a need for an alternative method to the 
Hilbert Transform in constructing an energy-frequency-time 
distribution. The Hilbert Transform is an intensive calcula- 
tion process that requires powerful computers beyond many 
end users, including ordinary engineers and scientists, 55 
resources. Furthermore its results are not intuitive. If a signal 
to be analyzed has a surge-like behavior, i.e. the signal 
contains a very high amplitude for a short period time such 
as sound signal of a gunshot, the method disclosed in the 
above references may not generate reasonable results due to 60 
limitations imposed by spline fitting, which is used in the 
EMD in constructing extrema envelopes of the signal. The 
problem with spline fitting applied to the situation just 
mentioned, might result in both overshoot/undershoot prob- 
lems and eventually divergence of IMF’S. 65 

Zero Crossing points are the point where the voltage 
polarity of a waveform changes from negative to positive (or 


where the expected number of extrema per unit time is given 
by 
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in which m i is the ith moment of the spectrum Therefore, the 
parameter, v, defined as 
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offers a standard bandwidth measure. For a narrow band 
signal v=0, the expected numbers of extrema and zero 
crossings have to equal. However, the previous equation 
defines the bandwidth in the global sense; it is overly 
restrictive and lacks precision at the same time as there 
might be more than two extrema between two zero crossing 
points. Consequently, the bandwidth limitation on the Hil- 
bert transform to give a meaningful instantaneous frequency 
has never been firmly established. 

SUMMARY OF THE INVENTION 

The present invention is directed to a computer imple- 
mented method of analyzing a physical signal. The method 
includes the steps of inputting the signal, extracting a set of 
Intrinsic Mode Functions from the signal, generating a set of 
mean frequency functions from the Intrinsic Mode Func- 
tions, and generating the instantaneous frequency based on 
critical points of the signal. 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIGS. l(a)-(£>) are high-level flowcharts describing the 
overall inventive method of extracting Intrinsic Mode Func- 
tions (IMFs) via Empirical Mode Decomposition (EMD) 
and using Generalized Zero Crossing to obtain the Instan- 
taneous Frequency. 

FIG. 2 is a graph illustrating an upper envelope, a lower 
envelope, an envelope mean and an original signal data. 

FIG. 3(a) is a graph of a sound signal of a gunshot. 

FIG. 3(b) is a graph describing the steps of constructing 
upper and lower envelopes of a part of the signal of FIG. 
3(a) with cubic spline fitting and generating an envelope 
mean from the envelopes. 

FIG. 3(c) is a graph describing the steps of constructing 
upper and lower envelopes of a part of the signal of FIG. 
3(a) with straight line fitting and generating an envelope 
mean from the envelopes. 
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FIG. 3(d) is a graph displaying a collection of Intrinsic 
Mode Functions (IMFs) extracted via Empirical Mode 
Decomposition (EMD) with cubic spline fitting as in FIG. 

m 

FIG. 3(e) is a graph displaying a collection of Intrinsic 
Mode Functions (IMF) extracted via Empirical Mode 
Decomposition (EMD) with straight line fitting as in FIG. 
3(c). 

FIGS. 4 (a)-(b) are graphs of data describing the Gener- 
alized Zero -Crossing approach. 

DETAILED DESCRIPTION OF PREFERRED 
EMBODIMENT OF THE INVENTION 

FIG. 1 (a) illustrates the overall inventive method of 
extracting IMFs via EMD including the Sifting Process in 
steps 220 through 270. First, the signal can be inputted in 
step 100. 

After inputting the signal in step 100, the analog signal is 
converted to the digital domain suitable for computer pro- 
cessing in the A/D conversion step 105. Depending upon 
whether the input signal is analog or digital step 105 may be 
bypassed. Thereafter, the Sifting Process (steps 107 through 
270) is applied to Sift the signal with the Empirical Mode 
Decomposition method and thereby extract the intrinsic 
mode function(s). 

The signal from step 100 is first windowed by framing the 
end points in step 107. Then, the Sifting Process begins at 
step 200 by identifying local maximum values of the digi- 
tized, framed physical signal from step 107. 

Then, the method constructs an upper envelope 20 of the 
signal 10 in step 210. The upper envelope 20 is shown in 
FIG. 2 using a dot-dashed line. This upper envelope 20 is 
constructed with a cubic spline that is fitted to the local 
maxima. 

Next the local minimum values of the signal 10 are 
identified in step 220. To complete the envelope, a lower 
envelope 30 is constructed from the local minimum values 
in step 230. The lower envelope 30 is also shown in FIG. 2 
using a dot-dash line. Like the upper envelope 20, the lower 
envelope 30 is constructed with a cubic spline that is fitted 
to the local minima. 

The upper and lower envelopes 20, 30 should encompass 
all the data within the physical signal 10. From the upper and 
lower envelopes 20, 30, an envelope mean 40 is the deter- 
mined in step 240. The envelope mean 40 is the mean value 
of the upper and lower envelopes 20, 30. As can be seen in 
FIG. 2, the envelope mean 40 bisects the physical signal 10 
quite well. 

Referring to steps 300 through 330, Generalized Zero- 
Crossing (GZC) is a method of computing the Instantaneous 
Frequency for a nonlinear and nonstationary signal based on 
full and partial periods around a data point 8 of the signal as 
shown in FIG. 4. Zero crossing points are points at which a 
signal under consideration crosses a zero value. 

In the GZC approach, all zero -crossings 2 and local 
extrema 4, 6 are defined as critical points. As shown in the 
FIG. 4 (a), the local extrema consist of local maxima 4 and 
local minima 6. For this approach, the signals to be analyzed 
are only meaningful Intrinsic Mode Functions (IMFs) when 
obtained via Empirical Mode Decomposition (EMD). 

The time intervals between all the combinations of the 
critical points 2, 4, 6 are considered as the whole or partial 
wave period (T/, T 2 y , and T 4 ). For example, the period 
between two consecutive up (or down) zero -crossings or two 
consecutive maxima 4 (or minima 6) can be counted as one 
period. Each given point along the time axis will have four 
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different values from this class of period, designated as T 1 . 
Then, the period between consecutive zero -crossings (from 
up to the next down zero -crossing, or from down to the next 
up zero -crossing), or consecutive extrema (from maximum 
5 4 to the next minimum 6, or from minimum 6 to the next 
maximum 4) can be counted as a half period. Each point 
along the time axis will have two different values from this 
class of period, designated as T 2 . Finally, the period between 
one kind of extrema to the next zero-crossing point, or from 
10 one kind of zero -crossings to the next extremum can be 
counted as a quarter period. Each point along the time axis 
will have one value from this class of period designated as 
T 4 . Clearly, the quarter class, T 4 , is the most local and the 
full wave class, T 1? is the least local. In total, there exist 
15 seven different period values at each point along the time 
axis. As T 4 , T 2 , and T ± are weighted by factors of 4, 2 and 
1, respectively, the weighted mean frequency at each point 
along the time axis can be computed as 

20 

_ 1 / 1 ( 1 1 W 1 1 1 1 Yl (1) 

w l2\T4 + l^ 2 1 + T 2 2 J + l^ 1 1 ' f ^ 1 24 '^ 1 3 + 7?JJ , 


25 wherein 

co is mean frequency; 

T ± x are full periods (x=l, 2, 3, and 4) enclosing the point 
under consideration; 

T 2 y are half periods (y=l and 2) enclosing the point under 
30 consideration; and 

T 4 is a quarter period enclosing the point under consid- 
eration; 

An alternative is to compute a non-weighted as 

35 

_1M_ (j_ J_ j_ J_)\ ( 2 ) 

m ~ 7\4 T 4 + U T\ + 2T}) + {t} + T\ + T\ + 7f )) 

40 It is to be noted that various combinations of the critical 
points and weights can be used to compute the mean 
frequency. 

Referring to step 300, a mean frequency function 42 of the 
signal 41 can be obtained based on a collection of the mean 
45 frequencies along the time axis as shown in FIG. 4. This 
process continues until all of the mean frequency functions 
are obtained for all of the IMFs as in steps 300 through 305. 
Summing up the mean frequency functions (step 310) gen- 
erates a function of frequency versus time, which is the 
50 instantaneous frequency. 

This approach is the most direct, local, and also the most 
accurate in the mean. Furthermore, this approach will also 
give a statistical measure of the scattering of the frequency 
value. 

55 This approach provides its crude localization, only down 
to a quarter wavelength. This approach tends to display its 
inability to represent the detailed waveform distortion; it 
admits no harmonics. Unless the waveform contains asym- 
metries (either up and down, or left and right), the GZC will 
60 give it the same frequency as a sinusoidal wave. For most 
practical applications, this mean frequency localized down 
to quarter a wave period is already a well-accepted result. As 
this method physically measures the period, or part thereof, 
the values obtained can serve as the best local mean over the 
65 period to which it applies. 

Referring to step 350 of FIG. 1(a) , serious problems of the 
spline fitting can occur near the ends, where the cubic spline 
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being fitted can have large swings. The end swings, if not 
fixed, can eventually propagate inward and corrupt the 
whole data span especially in the low frequency compo- 
nents. For example, if the signal contains a very high 
amplitude for a short period time such as the sound signal of 5 
a gunshot as shown in FIG. 3(a), the method disclosed in the 
above references may not generate reasonable results due to 
extreme overshoot in the cubic spline fitting. Referring to 
FIG. 3(c), which deals with a problematic area of the data of 
FIG. 3(a), the upper and lower envelopes 21, 31 are con- 10 
structed by connecting local maxima and minima with 
straight lines, respectively. For comparison, FIG. 3(b) shows 
the upper and lower envelopes 21', 31', which are con- 
structed with the cubic spline fitting. As anticipated, the 
EMD based on the cubic spline fitting does not generate 15 
reasonable IMFs as shown in FIG. 3(d) whereas IMFs 
obtained through the straight line fitting gradually converge 
to a monotonically decreasing signal as shown in FIG. 3(c). 

It is to be understood that the above -described arrange- 
ments are only illustrative of the application of the principles 20 
of the present invention. Numerous modifications and alter- 
native arrangements may be devised by those skilled in the 
art without departing from the spirit and scope of the present 
invention and the appended claims are intended to cover 
such modifications and arrangements. 25 

What is claimed is: 

1. The computer implemented method of analyzing a 
physical signal from a physical device comprising the steps 
of; 

a. inputting the physical signal; 30 

b. extracting a set of Intrinsic Mode Functions from the 
physical signal; 

c. generating a set of mean frequency functions from the 
Intrinsic Mode Functions, wherein the step of generat- 
ing a set of mean frequency functions includes com- 35 
puting the mean frequency at a point along the time 
scale for one of the Intrinsic Mode Functions and 
continuing to perform the computing step for all of the 
Intrinsic Mode Functions; and, 

d. displaying said set of mean frequency functions. 40 

2. The computer implemented method as in claim 1, 
wherein the mean frequency at a point under consideration 
is a weighted mean frequency. 

3. The computer implemented method as in claim 1, 
wherein the extracting a set of Intrinsic Mode Functions 45 
from the physical signal comprises: 

recursively sifting the physical signal via Empirical Mode 
Decomposition to extract an intrinsic mode function 
indicative of an intrinsic oscillatory mode; 
generating a residual signal by subtracting the intrinsic 50 
mode function from the physical signal; 
treating the residual signal as the physical signal during a 
next iteration of said recursive sifting step; and 
iterating to perform said recursive sifting to generate an 
n-th intrinsic mode function of an n-th intrinsic oscil- 55 
latory mode until a stopping condition is met. 

4. The computer implemented method of analyzing a 
physical signal according to claim 3, wherein said recursive 
sifting includes: 

identifying local maximum values in the physical signal; 60 
constructing an upper envelope of said physical signal 
from the identified local maximum values; 
identifying local minimum values in said physical signal; 
constructing a lower envelope of said physical signal from 
identified local minimum values; 65 

determining an envelope mean from the upper and lower 
envelopes; 
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generating a component signal by subtracting the enve- 
lope mean from said physical signal; 
treating the component signal as the physical signal; and 
recursively performing said sifting until successive com- 
ponent signals are substantially equal. 

5. The computer implemented method of analyzing a 
physical signal according to claim 4, wherein the step of 
constructing a lower envelope of the physical signal includes 
connecting the identified local minimum values with straight 
lines; and the step of constructing an upper envelope of the 
physical signal includes connecting the identified local 
maximum values with straight lines. 

6. The computer implemented method of analyzing a 
physical signal according to claim 4, wherein the step of 
constructing a lower envelope of the physical signal includes 
connecting the identified local minimum values with cubic 
spline fitting; and the step of constructing a upper envelope 
of said physical signal includes connecting the identified 
local maximum values with cubic spline fitting. 

7. The computer implemented method as in claim 1 
further comprising: the step of summing up the mean 
frequency functions. 

8. The computer implemented method as in claim 7 
further comprising the step of: 

displaying the sum of the mean frequency functions. 

9. A computer implemented method of analyzing a physi- 
cal signal from a physical device comprising the steps of: 

a. inputting the physical signal; 

b. extracting a set of Intrinsic Mode Functions from the 
physical signal; 

c. generating an instantaneous frequency based on critical 
points of the signal by generating a set of mean 
frequency functions from the Intrinsic Mode Functions, 
wherein the step of generating a set of mean frequency 
functions includes computing the mean frequency at a 
point along the time scale for one of the Intrinsic Mode 
Functions; 

d. continuing to perform the computing step for all of the 
Intrinsic Mode Functions; and, 

e. displaying said instantaneous frequency. 

10. The computer implemented method as in claim 9, 
wherein the mean frequency at a point under consideration 
is a weighted mean frequency. 

11. The computer implemented method as in claim 9, 
wherein extracting a set of Intrinsic Mode Functions from 
the physical signal comprises: 

recursively sifting the physical signal via Empirical Mode 
Decomposition to extract an intrinsic mode function 
indicative of an intrinsic oscillatory mode; 
generating a residual signal by subtracting the intrinsic 
mode function from the physical signal; 
treating the residual signal as the physical signal during a 
next iteration of said recursive sifting step; and 
iterating to perform said recursive sifting to generate an 
n-th intrinsic mode function indicative of an n-th intrin- 
sic oscillatory mode until a stopping condition is met. 

12. The computer implemented method of analyzing a 
physical signal according to claim 11 , wherein said recursive 
sifting including: 

identifying local maximum values in the physical signal; 
constructing an upper envelope of the signal from the 
identified local maximum values; 
identifying local minimum values in the physical signal; 
constructing a lower envelope of said physical signal from 
the identified local minimum values; 
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determining an envelope mean from the upper and lower 
envelopes; 

generating a component signal by subtracting the enve- 
lope mean from said physical signal; 
treating the component signal as the physical signal; and 5 
recursively performing said sifting until successive com- 
ponent signals are substantially equal. 

13 . The computer implemented method or analyzing a 
physical signal according to claim 12, wherein the step of 
constructing a lower envelope of the physical signal includes 10 
connecting the identified local minimum values with straight 
lines; and the step of constructing an upper envelope of the 
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physical signal includes connecting the identified local 
maximum values with straight lines. 

14 . The computer implemented method of analyzing a 
physical signal according to claim 12, wherein the step of 
constructing a lower envelope of the physical signal includes 
connecting the identified local minimum values with cubic 
spline fitting; and the step of constructing a upper envelope 
of said physical signal includes connecting the identified 
local maximum values with cubic spline fitting. 



